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INTRODUCTION 

■  '  ■••••  v  i ion 

The  moving  finite  element  (MFE)  method  is  a  new  approach  for  numerically 
solving  partial  differential  equation  (PDE)  systems; U»2»3)  is  particularly 
well  suited  for  resolving  PDE  solutions  which  may  contain  large,  multiple 
gradients  over  highly  disparate  scales  in  both  space  and  time.  These  types 
of  POE's  abound  in  such  basic  technical  disciplines  as  aerodynamics  (with  em¬ 
phasis  on  shear  layers,  shocks  and  their  possible  interactions),  combustion, 
plasma  physics,  material  interface  phenomena,  continuum  mechanics,  and  other 
transport  processes. 

In  the  MFE  method,  grid  co-ordinates  themselves  are  dependent  variables 
which  are  calculated  continuously  at  each  time  step  in  order  to  minimize  PDE 
residuals.  This  feature  has  successfully  suppressed  numerical  dissipation  to 
very  low  levels  and  has  resolved  accurately  in  1-D  those  physical  transport 
processes  which  may  occur  over  extremely  small  scales  simultaneously  with 
other  physical  processes  which  may  span  macroscopic  scales.  The  objectives 
of  the  presently  reported  research  effort  were  to  extend  the  MFE  method  to 
2-D  and  to  investigate  its  basic  properties  and  needs  for  solving  important 
PDE's  in  2-D.  This,  of  course,  is  an  on-going  effort  which  is  in  only  the 
early  stages  of  discovery  and  development.  This  report  describes  the  work 
tasks  and  results  which  have  contributed  to  attainment  of  the  objectives 
cited  above.  The  results  of  work  performed  to  date  can  be  summarized  briefly 
as  an  introduction  to  the  more  detailed  reporting  which  follows  in  subsequent 
sections. 

It  was  established  during  the  first  year  of  this  study  that  the  MFE 
method  does  indeed  extend  logically  and  practically  to  2-D.  This  fact  has 
been  substantiated  in  recent  work  by  others  as  well.5  An  initial  experi¬ 
mental  2-D  MFE  code  version  was  developed  in  the  first  year  of  our  work 
under  AFOSR  support  for  the  purpose  of  conducting  continuing  theoretical 
and  applied  mathematical  research  In  diverse  scientific  contexts.  Also  in 
the  first  year  of  study,  effective  MFE  node  movement  properties  and  signi¬ 
ficant  node  savings  were  demonstrated  for  simple  —  but  yet  significant  — 


problems  In  2-0.  The  MFE  code  structure  was  found  to  be  amenable  to  vector- 
izatlon  and  use  on  envisioned  advanced  computers.  Work  during  the  second 
year  of  study  focussed  on  such  essential  issues  and  research  needs  as  ODE 
Integrators  for  numerical  POE  solution  methods,  MFE  regularization  functions, 
grid  node  biasing  effects,  and  linear  solvers  for  large  skewed  matrices  which 
do  not  contain  diagonally  dominant  terms.  Work  in  the  third  year  has  ad¬ 
dressed  the  possible  resolution  of  highly  disparate  PDE  scales  in  order  to 
analyze  physical  dissipation  effects  in  shocks,  distinct  from  artificial  or 
purely  numerical  diffusion  effects.  This  has  involved  development  in  the  MFE 
method  of  physical,  non-slip  boundary  conditions  and  of  improved  numerical 
conditioning  of  matrix  solvers  for  the  discretized  equations  of  the  MFE 
method.  Detailed  accounts  of  research  during  this  three  year  project  are 
presented  in  the  report  which  follows. 


Mathematical  Background  of  the  MFE  Method 


In  order  to  present  a  coherent  report  of  work  to  the  broadest  range 
of  readers,  we  present  here  a  brief  mathematical  sketch  of  the  MFE  method. 
Consider  a  general  system  of  partial  differential  equations  (PDE's),  U  =  L(U) 


“1  *  «-x(UJ  .  (lJ 

"N  *  MU)  . 

Using  piecewise  linear  approxminants  of  Uj  ...  u^,*  which  are  of  the 
form  u  *  mx  +  ny  +  p,  on  a  hexagonal ly  connected  triangular  mesh  IFigures 
1-4),  application  of  the  chain  rule  to  the  differentiation  of  u^  gives: 

u.  =  £ak.akJ  +  x-8kJ  +  y.yk’5  ,  where  U) 


3u.  3u,  3u. 

_ L.  .  ftkJ  * _ •  vifJ  - _ - 

aakj  •  Bk  3x4  •  Yk  ayj 


The  key  aspect  of  the  MFE  method  lies  in  the  rigorous  retention  of  the 

•  • 

nodal  motion  terms  Involving  x  and  y.  In  most  conventional  PDE  solution 
methods,  the  nodal  motion  terms  are  either  set  to  zero  (fixed  nodes)  or 
allowed  to  assume  some  arbitrary  selected  values  (e.g.,  mean  fluid  velocities 
in  Lagrangian  codes). 

Alternatively,  the  MFE  method  uses  a  very  fundamental  criterion  for  the 

evaluation  of  grid  node  motions  (x  and  y)  at  every  time  step.  That  is,  the 

basic  MFE  method  is  formulated  by  requiring  that  all  of  the  time  derivatives 

.  ... 

aIj,  ....,  aNj,  xj,  yj  be  found  at  each  instant  in  such  a  way  as  to  minimize 
the  L2  norm  of  the  PDE  residual,  ||  U  -  L(U)  |j  .  This  has  been  found  to 
minimize  PDE  solution  errors  and  to  dramatically  reduce  the  number  of  grid 


♦The  present  mesh  triangulation  has  been  chosen  for  simplicity.  Many  other 
grid  meshing  schemes  are,  of  course,  possible;  and  many  of  the  alternative 
schemes  should  be  investigated  in  future  work.  It  is  interesting  to  note, 
however,  that  numerous  MFE  results  confirm  that  PDE  solution  approximants 
of  higher  degree  are  not  nearly  so  critical  for  attaining  high  levels  of 
stability,  accuracy  and  resolution  in  an  optimal  moving  node  solution  method 
as  they  are  in  fixed  node  or  less  optimal  adaptive  PDE  solution  methods. 


I 


Figure  4.  MFE  Nodal  Grid. 


\l 


nodes.  In  order  to  handle  the  degeneracies  which  can  sometimes  occur  in  such 
a  parametrization  of  PDE  solutions,  regularization  techniques  are  applied. 
That  is,  the  more  general  function. 


*  =  »1  II  “J  -  II  2  +  •••  “n  1 1  “H  -  ln(u) 


2  *  triangle  '  s/ 

alti tudes 


is  minimized  with  respect  to  the  parameter  derivatives  al . aN.,  x.,  y. 

J  J  J  J 

The  present  code  applies  the  most  elementary  form  of  regularization  functions 
in  which  the  expressions  for  e  -2 (d  - )  and  e.S  .  *  e  -  (d  .  )S  .(d.)  become  extremely 

J  tl  J  J  J  J  J  J 

large  and  thus  restrict  node  motions  when  any  triangle  altitude  dj  approaches 
a  user-specified  minimum  separation  tolerance.  Regularization  techniques 
also  serve  to  move  nodes  according  to  certain  physical  criteria  in  addition 
to  providing  numerically  well-conditioned  node  movement  properties.  The 
variational  equations  for  this  minimization  procedure  are  then  represented 
by  the  large  set  of  OUt's  which  follows: 


5Z[(cr*,  ak^)akj+  (&c^,  +  (yk^,  a’jy^J  * 

J 


(Lk(U),  a* )  for  k  *  1,  ...  N, 


-5- 


(4b) 


Z)  w2  Bk^ak.  +  (Bkj,  Bk^Jx.  +  (Ykj»  Bk1  )y  .  1  « 

k»l  k  j  L  J  J 


No  j 

Y*  w  (L. (U),  Bk  )  +  (regularization  terms) 
k=l  k  K 


Z)  w2  Z3  [(<***,  Yk^ak.  +  (Bk^,  Yk*)x.  +  (Yk**,  Yk1)^.]- 
k=l  k  j  L  J 

(4c) 

N 

13  w2  (L.  (U),  Yk1)  +  (regularization  terms) 
k*l  k  K 

The  sums  on  j  in  Equations  (4)  run  over  the  seven  neighboring  nodes  of 
i  (including  the  ith  node  itself)  in  the  hexagonal  grid.  Equations  (4)  thus 
provide  the  basic  working  equations  of  the  MFE  method  in  2-1).  This  system  of 


OUE's  is  of  the  form 

R(y)  =  C(y)y  -  g(y)  *  0  ,  (5) 

where  y(t)  *  (alj,  . xi,  yi;  al2»  ...»  X2.  y2i  . )  is  the  vector 


of  unknown  parameters,  and  the  "mass  matrix"  C(y)  is  symmetric  and  positive 
definite.  This  system  of  ODE's  can  be  quite  stiff,  and  stiff  ODE  solution 
methods  are  thus  required  in  order  to  effectively  obtain  the  numerical  solu¬ 
tions  of  Equations  (4)  and  (5).  Research  which  attempts  to  eliminate  this 
stiffness  feature  is  In  a  very  early  stage  of  investigation, and  preliminary 
results  are  not  reported  at  this  time. 

The  stiffly  stable  numerical  integration  methods  which  can  be  used  to 
solve  the  system  of  ordinary  differential  equations  (4)  and  (5)  usually  con¬ 
tain  a  series  of  backward  Cauchy-Euler  steps  interspersed  with  interpolations 
and  extrapolations,  all  amidst  error-controlling  tests  and  strategies,  a 
backward  Cauchy-Euler  (BCE)  step  is  obtained  by  replacing  y  in  Equation  (5) 
by  the  backward  difference  (y  -  y)/At,  where  y  is  the  known  parameter  vector 
at  the  previous  time,  and  y  Is  the  unknown  parameter  vector  at  the  current 
time.  Upon  linearization,  the  BCE  equation  then  reads  as 


R(y,  y)  s  R(y)  +  R  (y)  •  6y  *  o  .  lb) 

The  structure  of  the  mass  matrix  C(y)  determines  in  large  measure  the 
degree  of  difficulty  and  computational  cost  of  obtaining  numerical  solu¬ 
tions  for  6y  in  Equation  (6)  and  thus  for  the  y  array  in  Equations  (4)  and 
(5).  Both  an  implicit  Runge-Kutta  method  (DIRK2)  and  the  Gear  method  are 
used  interchangeably  in  the  present  test  MFE  code  in  2-0;  both  of  these 
stiff  ODE  solvers  incorporate  the  backward  Cauchy-Euler  steps  described 
above. 

The  structure  of  the  mass  matrix  C(y)  can  be  seen  by  considering  the 
15th  node  in  a  6  x  6  grid  mesh  of  the  type  shown  in  Figures  1-4.  The 
coupling  among  nodes  in  Equation  (4)  is  represented  schematically  by 

20  21 

V*  _ 15  SSS\  16 


and  the  block  structure  of  the  mass  matrix  is  shown  in  Figure  5.  For  a  PUE 
system  with  2  "problem  variables,"  the  segment  of  the  dependent  variable 
array  y  associated  with  the  15th  node  is  (y)15  *  lal15,  a215»  xis»  y15)*» 
and  the  4  x  4  block  C15  16  is 


(°16’  °15) 


‘16’  15' 


(6l16’  a15J 

(rli6.  «15 

(6216’  “IS* 

(y^16’  ”15 

15,16  ^al6‘  6l15^  (°16’  ®215^  6I<15^  ^/Ykl6’  6I<15^ 

2  2 

(a16*  yl15*  {o16’  y215*  Yk15*  Yk15^ 


♦The  co-ordinate  variable  yj5  of  the  15th  MFE  node  is  not  to  be  confused 
with  the  dependent  variable  of  the  MFE  method  which  is  also  denoted  by  y. 


Although  direct  L-U  decomposition  and  solution  of  Equation  (6)  Is  rela¬ 
tively  slow.  It  Is  reliable  and  was  thus  used  In  the  Initial  code  versions  In 
the  presently  reported  work.  More  efficient  matrix  solution  methods  remain 
under  investigation. 


•••  9  10  11  12  13  14  15  16  17  18  19  20  21 
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Figure  5.  Block  Structure  of  the  Mass  Matrix  C(y). 

A  large  amount  of  the  computational  effort  in  the  MFE  method  is  devoted 
to  repetitive  evaluations  of  the  inner  product  terms  in  Equations  (4),  which 
are  then  used  iteratively  in  the  numerical  ODE  integration  procedures,  be¬ 
cause  invariant  geometrical  relationships  can  be  exploited  frequently  in  the 
evaluation  of  inner  products,  these  geometrical  relationships  are  generated 
and  encoded  once,  and  for  all  time,  from  the  initialization  data  at  the 
outset  of  code  calculations.  Figures  6-8  illustrate  the  grid  mesh  conventions 
which  would  be  established  and  encoded  by  the  present  2-D  research  code  for 
42  nodes  on  a  7  x  6  grid.  The  manner  in  which  these  relationships  are  used 
can  be  seen  readily  by  considering  a  single  conservation  equation 

Ut  *  -fx  “  9y  +  Uxx  +  Uyy  (b) 

and  the  piecewise  linear  basis  functions  about  the  ith  no<je  0f  the  form 


aj  *  ajx  +  bjy  +  cj 
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Figure  6.  Triangles. 
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The  inner  products  (aj,  ctj)  which  appear  in  the  mass  matrix  of  Equations 
(4)  are  evaluated  as 


a1 


.  «,)  -  / 


6T 


dxdy  (a^x  +  b.y  ♦  c^)' 


(10) 


where  the  ith  node  is  surrounded  by  six  adjoining  triangles  T.  The  two- 
dimensional  integration  over  a  triangle  T  in  Equation  (10)  can  be  performed 
either  analytically  or  by  the  midpoint  rule  according  to  the  formula 


Jw(x,  y)  dxdy  z  (A/3)  [w(P1)  +  w(P2)  ♦  w(P3)]  , 


(11) 


where  A  is  the  area  of  triangle  T,  and  P^  p2,  and  P3  are  the  midpoints  of 
the  sides.  This  simple  midpoint  formula  of  integration  is  exact  whenever  the 
integrand,  w(x,  y)  is  a  quadratic  function  in  x  and  y.  The  present  research 
code  uses  the  above  midpoint  rule  and  contains  options  for  the  use  of  a  com¬ 
posite  midpoint  and  some  analytic  integration  methods,  as  well.  Using  the 
relationships  between  basis  functions  and  Yj  *  UyCtj,  and  noting 

that  ux  and  uy  are  constants,  the  inner  products  (aj,  8i),  (81.  8^),  (8*,  yj) 
and  (y j ,  y.|)  can  also  be  evaluated  by  the  midpoint  rules.  For  a  general 
operator  L(u)  =  f,  it  can  be  shown  that 

("fx,  cij)  =I](ax)1  dxdy  ft  and  {Vi) 

|-fx.  8,)  ‘£»x  |-(<>x>,  fixtof  ♦  ^  »jjf  (13) 

edges  of  T 

The  integration  on  ds  is  performed  with  respect  to  arc  length  s  along  a 
triangle  edge  and  *  t(s)  is  the  linear  function  of  s  which  rises  from  a 
value  0.  at  an  outer  node  to  a  value  1.0  at  the  ith  central  node.  Edge 
integrals  can  be  evaluated  analytically  or  by  numerical  quadrature.  Simpson's 
rule  and  a  composite  Simpson's  rule  are  available  for  use  in  the  present 
research  code.  The  derivatives  ax,  ux,  and  the  x  component  of  the  unit 
outward  normal  nj  are  readily  determined  at  all  times  by  Invariant  formulae 
from  the  known  amplitudes  and  nodal  co-ordinates  at  triangle  vertices. 


Inner  products  of  Laplacian  operators  are  more  complicated  to  evaluate, 
requiring,  in  addition  to  the  quantities  above,  the  unit  tangent  vectors  and 
their  x  and  y  components. The  formulae  for  these  quantities  are  also 
encoded  at  the  outset  of  computations  for  repeated  use  in  later  code  calcu¬ 
lations. 


Figure  9 


Schematic  Representation  of  Major  Functions  Performed  in  the 
Present  MFE  Codes  in  Both  1-D  and  2-0. 


The  overall  2-1)  code  structure  is  conceptually  simple,  as  Indicated  In 
Figure  9.  The  detailed  code  structure  contains  an  assembly  of  approximately 
forty  subroutines  which  perform  modularly  the  numerous  substasks  which  are 
required  In  order  to  execute  the  major  code  functions  Indicated  In  Figure  9. 
This  structure  Is  Intended  to  facilitate  mathematical  research  on  alternative 
(actually.  Interchangeable)  ODE  solution  methods,  matrix  solution  methods, 
node  control  strategies,  boundary  conditions,  and  such  other  tasks  as  grid 
mesh  generation  and  Issues  of  numerical  analysis  which  must  be  Investigated 
In  new  POE  solution  methods.  This  code  structure  of  the  MFE  method  also 
appears  to  be  highly  compatible  with  vector  and  parallel  processing  computers. 

Because  the  mass  matrix  C(y)  In  Equation  (5)  can  become  singular  (when¬ 
ever  all  components  u  have  a  flat  portion  in  their  graph  at  some  node  (x^.y^)), 
regularization  terms  are  included  in  the  minimization  problem  for  the  x  and  y 
equations.  The  minimization  problem  for  x  and  y  reads  =  0  and  =  0  , 
where  X|c  k 


-Lf  (U)||  *  £  (tya,)  dj  -Sjtdj))2 

triangle 

altitudes 


This  minimization  yields  regularization  terms  of  the  form: 


in  the  x  equation 


e ,/e  .d .  -  S  ,(d)\ 
J\  0  0  J  J 


(15a) 


in  the  y  equation 


fe.d.  -  S.(d.)).«i 

l  0  J  J  J  J  3y 


(15b) 


It  follows  that  the  terms  d 


3d  •  _  3d 

*  %  *f  —  • 

dx.  k  dy. 


y  and  the  terms  and  ~T 
k  3x.  3y 


In  Equations  (15)  are  dependent  primarily  upon  grid  geometry;  and  their  in¬ 
variant  formulae  are  thus  generated  at  the  outset  and  encoded  for  repeated 
calculations  later  in  the  code  in  a  manner  similar  to  many  portions  of  the 
residual  evaluations  which  were  discussed  above. 

The  functions  e1 2 3  and  eS  which  appear  in  Equations  (15)  are  treated  in 
exactly  the  same  manner  as  an  advanced  penalty  function  formulation  which 
has  been  tested  previously  in  1-0  investigations^*2^  of  effective  regular¬ 
ization  methods;  that  is. 


(d  -  SEPMIN) 


(d  -  SEPMIN; 

where  SEPMIN  is  a  minimum  triangle  altitude.  A  detailed  node  controlling 
policy  has  been  developed  for  the  selection  of  values  for  the  constants  CA 

and  C2  in  terms  of  local  truncation  errors  and  of  the  inner  products  ( 6, 6) 
and  (y,  y)  which  appear  in  the  *  and  $  equations.  A  discussion  of  this 
policy  in  specific  1-0  problem  applications  appears  in  References  1  ana  2. 

A  detailed  discussion  of  this  policy  for  2-0  applications  is  deferred  until 
several  additional  generalizations  are  discussed.  These  generalizations  will 
also  have  the  effect  of  preventing  automatically  the  inversions  (tangling) 
of  triangles.  Such  inversions  are  presently  detected  and  prevented  by 
testing  the  aspect  ratios  and  the  signs  of  triangle  determinant  quantities. 
Integration  time  steps  are  reduced  when  aspect  ratios  exceed  a  designated 
value  or  when  there  is  a  change  of  sign  in  the  determinant  quantities. 


Numerous  sample  problems  have  been  used  for  testing  of  the  MFE  method 
in  2-0.  These  problems  have  been  designed  to  test  such  numerical  aspects  of 
the  MFE  method  as: 


(1)  Inner  product  evaluations  (analytic  and  numerical  quadrature) 

(2)  ODE  integration  for  the  MFE  equations 

(3)  Matrix  solution  methods 
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(4)  Regularization  schemes  for  control  of  node  motions 

(5)  Jacobian  evaluations 

(6)  Boundary  conditions  (zero-Neumann  and  Dirichlet). 

The  POE's  in  these  test  problems  are  in  the  form  of  general 
conservation  equations, 

ut  =  "fx  -  9y  +  v(uxx  +  uyy}  * 

where  the  flux  functions  f  and  g  can  have  nearly  arbitrary  functional 
forms. 


In  addition  to  the  simple  examples  of  the  heat  equation  and  of  square 
wave  propagation  which  were  presented  in  Reference  3,  two  somewhat  more 
complex  test  examples  have  also  been  studied. 

Oblique  Propagation  of  a  Scalar  Wave 

This  example  considers  the  propagation  by  pure  advection  of  a  scalar 
wave  diagonally  across  the  initial  grid  mesh,  according  to  the  equation  ut  + 
ux  +  uy  *  0.  The  initial  conditions  for  this  example  are  shown  schematically 
in  Figure  10  and  are  expressed  as: 


u(x,y,0)  «  1.0 
u(x,y,0)  *  0. 
u(x,y,0)  ■  linear 


0.1  <  x  <  0.2;  0.1<y  <0.2 

x  <_  0.05;  x  0.25  and  all  y 
otherwise  . 


Dirichlet  conditions,  u(x,y,t)  =  0.,  are  applied  at  all  boundaries. 

The  co-ordinates  x  and  y  obey  the  following  Dirichlet  conditions:  x  =  0. 
along  the  y  axis;  x  ®  1.0  along  the  boundary  x  *  1.0,  all  y,  y  *  0.  along 
the  x  axis;  and  y  *  1.0  along  the  boundary  y  =  1.0,  all  x.  The  co-ordinate 
variables  obey  zero-Neumann  conditions  on  y  along  the  y  axis  and  along  the 
boundary  at  x  *  1.0  for  all  y  and  on  x  along  the  x  axis  and  along  the  bound¬ 
ary  at  y  «  1.0  for  all  x.  (That  is,  the  y  co-ordinate  is  free  to  slide  along 
the  vertical  boundaries,  and  the  x  co-ordinate  is  free  to  slide  along  the 
horizontal  boundaries.) 
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Using  a  6x6  grid  of  moving  nodes,  this  problem  is  run  from  t  *  0.  to 
t  ■  0.8,  which  spans  the  interval  of  free  propagation.  The  accuracy  of  the 
wave  profile  is  maintained  to  within  1  part  in  i03,  consistent  with  the 
local  truncation  error  constraint  in  the  ODE  solver.  The  velocity  of  pro¬ 
pagation  is  accurate  to  4  significant  figures.  The  mesh  is  observed  to  move 
flexibly  in  order  to  maintain  these  accuracies  throughout  the  problem  evolu¬ 
tion.  As  this  wave  approaches  the  upper  right-hand  corner  of  the  domain, 
aspect  ratios  of  some  of  the  grid  mesh  triangles  approach  values  on  the  order 
of  approximately  10 2,  with  no  adverse  effects.  The  triangle  aspect  ratios 
can  be  made  to  assume  values  which  are  an  order  magnitude  larger  oy  imposing 
Dirichlet  conditions  on  the  x  and  y  coordinates  at  the  boundaries.  It  was 
also  found  that  the  MFE  method  can  solve  this  problem  with  equal  facility 
and  efficiency  for  much  larger,  finite  gradients  of  the  scalar  wave  front. 

This  problem  can  readily  be  modified  so  that  the  scalar  wave  trajectory 
follows  a  circular  path  according  to  the  pure  advection  equation, 

u.  *  cos(t)  *  u  +  sin(t)  u  ,  -a  <  x  <  a,  -a  <  y  <  a  ,  ( 19 J 

t  x  y 

where  the  dimensions  a  are  sufficiently  large  to  contain  the  circular  tra¬ 
jectory.  In  this  MFE  solution,  the  scalar  wave  follows  its  circular  trajec¬ 
tory  accurately  and  returns  precisely  to  its  initial  position,  consistent 
with  the  local  truncation  error  constraint  in  the  ODE  solver,  after  a  com¬ 
plete  revolution  (t  *  2ir),  using  20  time  step  cycles.  The  grid  mesh  again 
expands  and  contracts  smoothly  and  flexibly  in  maintaining  the  desired 
solution  accuracy  at  all  times. 

Burger-Like  Equations. 

A  2-0  analog  of  Burger’s  1-0  model  equation  is  given  by  the  equations 

if  “ "  Hr*- "  iy (uv)  +  vv2u  (*u) 

~  Tx  (uv)  "  5y  (v2)  +  w2v  *  Ul) 
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Figure  10.  Initial  Conditions  for  Oblique  Propagation  of  Scalar  Wave 


Direction  of 
Propagation 
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Figure  11.  Transient  Propagation  of  a  Scalar  Wave. 
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where  u  and  v  can  be  viewed  as  x  and  y  components  of  a  fluid  velocity, 
respectively.  In  order  to  maintain  a  close  analogy  to  the  1-b  burger's  model 
results  which  were  discussed  in  Reference  2,  initial  and  boundary  conditions 
for  this  system  of  PDE's  are  first  selected  so  as  to  lead  to  the  propagation 
of  a  uniform,  step-like  wave  in  a  direction  parallel  to  the  x-axis;  that  is, 

u(x,y,Q)  =  0.  0.  £  x  £  1.0;  0.  £  y  £  1.0 

v(x,y,0)  *  1.0  0  £  y  1  0.100;  0  £  x  £  1.0 

v(x,y,0)  ■  0.  0.101  <  y  <  1.0;  0  <  x  <  1.0 

v(x,y,0)  *  linear  otherwise 

This  problem  is  solved  from  t  *  0  to  t  ■  0.5  on  a  4  x  19  grid  with  v»  10'*.* 
The  dependent  variable  v  obeys  zero-Neumann  conditions  along  the  y  axis  and 
along  the  (vertical)  boundary  at  x  3  1.0;  and  v  obeys  the  Oirichlet  condi¬ 
tions,  v  *  1.0  along  the  x  axis  and  v  *  0.  along  the  (horizontal)  boundary 
at  y  *  1.0.  All  interior  nodal  co-ordinates  obey  the  same  type  of  sliding 
boundary  conditions  as  were  used  in  the  scalar  wave  example  discussed  above. 
Figures  12-15  present  the  MFE  solutions  of  this  evolving  wave  front.  The 
extensive  migration  of  the  MFE  nodes  from  their  initial  positions  to  those 
positions  which  resolve  the  waveform  at  t*0.2  is  clearly  evident  in  Figures 
12  and  13.  The  speed  of  propagation  and  the  shock-like  wavefront  solutions 
are  resolved  to  accuracies  of  three  significant  figures,  or  better,  consis¬ 
tent  with  the  local  truncation  error  constraint  of  the  OUL  solver.  The 
magnitudes  of  the  wavefront  gradients  are  approximately  100  in  this  example, 
and  MFE  solutions  of  this  problem  can  be  obtained  with  similar  facility  and 
efficiency  for  much  larger  gradients  (corresponding  to  smaller  values  of  v 
in  Equations  (20)  and  (21)).  Consistent  with  arlier  results  in  Reference  7 
for  simple  square  wave  propagation  by  pure  advection.  It  is  found  also  in 
this  Burger-like  example  that  wide  latitude  can  be  exercised  in  the  selection 
of  node  controlling  parameters  in  the  functions  e*  and  eS  of  Equations  (16) 

*This  problem  can  be  solved  with  equal  effectiveness  on  an  MFE  grid  of  4xiu 
nodes.  The  4x19  grid  simply  represents  the  Initial  attempt  on  tnis  problem. 
The  figures  12-15  below  have  rotated  the  co-ordinate  axes  in  the  x-y  plane 
by  90*  from  the  conventional  orientation  (£  horizontal  and  y  vertical)  in 
order  to  Improve  the  viewing  angle  for  the  plotted  results  Tn  3-D.  The  terms 
"horizontal"  and  "vertical"  refer  strictly  to  the  conventional  orientation 
of  the  x-y  plane  throughout  this  discussion. 


u 


Tt»e  •  0 
«  k  19  MFE  Grid 
v  -  0.01 
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Figure  13.  MFE  Solution  of  Burger-like  Equation  in  2-D. 


and  (17).  When  relatively  weak  Internodal  forces  are  used,  no  nodal  devi¬ 
ation  (or  bias)  In  the  x  direction  Is  observed.  (There  are  no  transverse 
forces  In  the  PDE's,  per  se.)  As  the  Internodal  forces  are  Increased  to 
sufficiently  large  values,  the  nodes  can  be  forced  by  the  regularization 
terms  to  migrate  toward  the  x-axIs  for  the  nodal  triangulation  shown  In 
Figure  12 — meanwhile  maintaining  the  accuracies  mentioned  above.  The 
direction  of  such  forced  nodal  migration  can  be  reversed  by  using  the 
opposite  type  of  grid  triangulation,  and  this  node  migration  can  also  be 
altered  by  the  use  of  different  penalty  functions.  Node  control  is  thus 
very  flexible  and  desired  accuracies  are  readily  maintained.  When  PDE's 
contain  non-trlvlal  x-dependencies  In  their  operators,  the  PUE's  themselves 
resume  their  dominant  role  in  governing  the  positioning  of  the  nodes,  as 
will  be  seen  In  a  skewed  waveform  example  below.  Figure  15  shows  the 
formation  of  the  boundary  layer  at  the  right-hand  boundary,  as  this  MFE 
solution  approaches  the  correct  asymptotic  solution. 

An  example  of  a  skewed,  propagating  wavefront  (shock)  can  be  formulated 
in  terms  of  the  Burger-like  equations, 

ut  a  0  (22) 

vt  «  -1/2  (v2^  +  vV2  ,  (23) 

on  the  unit  square.  An  8xu  grid  of  MFE  nodes  is  used,  and  the  initial  con¬ 
ditions  for  u  and  v  on  uniformly  spaced  grid  nodes  are: 

u(x,y,0)  =  0  .  all  x,y 

v(x,y,0)  *  1+7t,  1+6t,  1+5t,...,1  at  nodes  1,  2,  3,...,B 
along  1st  row  (x-axis) 

v(x,y,0)  =  -1,  -(1+t),  -(1+2t),...,-(1+7t)  at  nodes  57,  5b, 
59,..., 64  along  top  row. 

The  Initial  values  of  v(x,y,0)  along  a  given  vertical  line  are  obtained 
by  linear  Interpolation  at  interior  nodes.  The  parameter  t  is  assigned  a 
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constant  value  of  0.01,  and  the  value  of  vin  Equations  (22)-(23)  Is  assigned 
a  value  of  0.01  In  the  present  run.  As  shown  In  Figure  16,  these  Initial 
conditions  on  v  simply  map  a  plane  In  which  v  has  the  values  +1.07  at  the  co¬ 
ordinates  (0,0);  +1.U  at  (1,0);  -1.0  at  (0,1);  and  -1.U7  at  (1,1).  Uirichlet 
boundary  conditions  are  maintained  on  v(x,y,t)  and  on  x  and  y  along  the  hori¬ 
zontal  boundaries;  zero-Neumann  boundary  conditions  are  applied  to  v(x,y,t) 
and  to  y  on  the  vertical  boundaries;  and  Dlrlchlet  boundary  conditions  are 
maintained  on  x  along  the  vertical  boundaries.  This  skewing  of  the  Initially 
counter-directed  velocity  components  along  the  top  and  bottom  boundaries 
leads  to  the  evolution  of  non-uniform  wavefront  solutions  which  are  seen 
in  the  results  below.  In  the  early  stages  of  solution,  prior  to  t  ■  O.b,  a 
projection  of  the  MFE  solution  on  the  x-y  plane  shows  two  counter-directed, 
quasi -horizontal  wave  impulses  which  propagate  from  top  to  bottom  and  from 
bottom  to  top  at  speeds  of  approximately  U.5.  At  t  «  u.5,  a  shock  is  formed 
when  the  propagating  impulses  encounter  each  other  near  the  horizontal  cen¬ 
terline  of  the  x-y  domain.  Subsequently,  a  skewed,  shock-like  waveform  is 
generated  and  propagates  in  the  serpentine  manner  shown  in  Figures  16-19.  The 
relatively  large  aspect  ratios  seen  in  Figure  19  for  the  MFE  mesh  at  t  =  2D. 
were  created  deliberately  by  the  use  of  Dlrichlet  boundary  conditions  on  the 
x  co-ordinates  along  the  x  axis  and  along  the  parallel  boundary  line  at  y  = 

1.  The  MFE  node  migration  was  fluid  throughout  and  exhibited  no  grid-biasing 
effects.  This  problem  was  run  from  t  *  U.  to  t  =  2U.  in  approximately  125 
time-step  cycles.  The  gradients  of  the  fully  developed  shock  are  on  the  order 
of  100,  consistent  with  the  present  value  of  v  *  0.01.  As  above,  MFt  solu¬ 
tions  of  this  problem  on  an  8  x  8  grid  can  be  obtained  for  much  larger  gra¬ 
dients  (smaller  values  of  v)  with  essentially  the  same  levels  of  robustness 
and  efficiency  as  are  seen  In  the  present  example. 


From  these  early  2-U  results  it  was  apparent  that  the  hexagonal ly 
connected  triangular  mesh  and  perhaps  several  other  possible  triangulation 
schemes  are  quite  compatible  with  the  MFE  method.  The  additional  degrees  of 
geometrical  freedom  which  are  available  for  error  minimizing  node  motions  in 
2-D  have  been  found  to  have  a  beneficial  effect  on  the  numerical  integration 
efficacy  of  the  MFE  method  in  2-D,  vis  i  vis  the  more  highly  constrained 
nodal  motions  in  1-D  MFE  solutions. (7.8)  y^e  nodal  movement  properties 
observed  in  these  Initial  2-U  results  thus  suggest  some  likely  implications 
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Figure  18.  MFE  Solution  of  Burger-like  Equations  for  a  Skewed  Waveform  in  2-D. 


Figure  19.  Projection  on  the  x-y  Plane  of  the  MFE  Solution  of  Burger-like 
Equations  for  a  Skewed  Waveform  in  2-D.  time  =20.;  8x8  MFE 
urid;  v  *  0.01. 


on  the  eventual  role  of  mesh  generation  needs  and  data  structure  Issues  In 
the  MFE  method.  So  long  as  the  MFE  method  exhibits  this  robustness  in  its 
continuous  node  movement  properties,  the  use  of  grid  mesh  generation  routines 
would  be  reduced  primarily  to  problem  initializations. 

In  order  to  study  the  MFE  nodal  migration  properties  under  more  demand¬ 
ing  circumstances,  the  initial  conditions  for  the  Burger's  model  above  can  be 
extended  so  as  to  create  a  more  highly  skewed  waveform  in  two  directions.  The 
PDE's  for  this  skewed  Burger's  model  problem  are  given  by: 

ut  *  -uux  -  vUy  +  v{uxx  +  uyy)  (24) 

vt  =  -uvx  -  vvy  +  v(vxx  ♦  vyy )  (25) 

where  u  is  the  x-component  of  velocity  and  v  is  the  y-component,  and  v  is  an 
effective  diffusion  coefficient.  Shocks  are  generated  with  gradient  magni¬ 
tudes  on  the  order  of  1/v.  Initial  conditions  which  produce  a  doubly  skewed 
wavefront  profile  are  shown  schematically  in  Figures  20  and  21.  (The  counter- 
posed  initial  velocity  fields  are  designed  to  create  an  evolving  shock  profile 
which  is  skewed  in  both  the  x  and  y  components  of  velocity.)  Boundary  condi¬ 
tions  are  given  by: 

u(0,y)  =  u(l,y)  *  0 

vx(0.  y)  =  vx(l,  y)  *  0. 
u(x,  1)  =  0.2  sin  irx 

u(x,  0)  =  -.2  sin  irx 
v(x,  1)  =  -1.  +  0.2  cos  *x 
v(x,  0)  =  1.  +  0.2  cos  7TX 

The  MFE  nodes  are  fixed  by  zero  Dirichlet  conditions  along  the  top  and  bottom 
horizontal  edges  of  the  domain.  The  nodes  are  free  to  move  vertically  by 
symmetric  boundary  conditions  along  the  left  and  right  edges  of  the  domain. 

This  problem  can  be  solved  readily  by  perhaps  many  POE  solution  methods 
whenever  v  assumes  sufficiently  large  values.  For  example,  a  value  of  v  * 

0.02  produces  shock  gradients  on  the  order  of  102. 
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The  MFE  method  requires  only  an  8  x  8  grid  to  give  reasonably  accurate 
solutions  to  this  problem,  and  Figures  22  and  23  show  very  accurate  MFE  solu¬ 
tions  on  a  12  x  12  grid.  Here  Figure  22  presents  an  Isometric  view  of  the 
evolving  profile  of  the  y  component  of  velocity  at  t  *  3.0,  well  after  the 
shock  has  formed  and  after  the  wavefront  has  undergone  significant  shearing. 
The  x-component  of  velocity  Is  sufficiently  sheared  that  a  hidden  line  plot, 
which  is  not  yet  available,  is  required  for  easy  Interpretation  by  the  naked 
eye.  The  MFE  grid  nodes  have  migrated  extensively  from  their  initial  posi¬ 
tions  as  can  be  seen  in  Figure  23  which  represents  the  grid  mesh  projected 
onto  the  x-y  plane  at  t  *  3.0.  Figures  24  and  25  present  contour  plots  for 
selected  constant  values  of  u  and  v,  respectively,  at  t  *  3.0.  It  is  evident 
from  the  magnitudes  of  shock  gradients  and  from  the  regions  of  significant 
curvature  which  span  nearly  the  entire  domain  that  an  alternative  POE  method 
with  a  fixed  grid  may  require  on  the  order  of  104,  or  more,  grid  nodes  in 
order  to  achieve  comparable  degrees  of  accuracy  in  this  problem. 

This  same  basic  problem  can  now  be  made  to  correspond  to  a  much  more 
demanding  physical  problem  by  letting  v  *  0.0U2.  Figure  2b  shows  an  isometric 
view  of  the  MFE  solution  on  a  16  x  16  grid  for  this  case.  Shock  gradients 
are  now  generated  with  magnitudes  of  several  times  103.  Before  discussing 
these  MFE  results  in  detail,  some  general  observations  should  be  discussed: 

It  is  likely  that  most  existing  PDE  method  using  either  a  fixed  grid  or  a 
less  than  optimal  adaptive  grid  would  require  on  the  order  of  105-10b  grid 
nodes  to  solve  this  problem  with  comparable  accuracy.  (We  note  as  an  aside 
that  numerous  inviscid  solvers  which  are  under  development  do  not  apply  at 
all  to  this  type  of  advection-diffusion  problem  because  the  Laplacian  is  an 
essential  mathematical  operator  whose  effects  must  be  rigorously  resolved  in 
advection-diffusion  PUE's.  Because  inviscid  solvers  do  not  generally  solve 
PDE's  which  contain  Laplacians,  they  generate  shocks  with  gradient  shapes  and 
magnitudes  that  are  governed  exclusively  by  the  selected  gridding  and/or  by 
the  purely  numerical  dissipative  processes  in  the  inviscid  method,  per  se. 
Consequently,  inviscid  solvers  would  have  no  chance  of  resolving  correctly 
any  of  those  physical  dissipation  effects  which  are  usually  expressed  by 
Laplacian  operators  and  are  present  with  fundamental  physical  significance 
In  transport  theory,  hydrodynamics,  plasma  physics,  continuum  mechanics,  and 
many  other  disciplines  in  the  physical  sciences.  This  critical  discussion 
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Figure  23.  MFE  grid  projections  on  the  X-Y  plane  at  t  *  3.0  In 
the  2-D  Burger-like  example  on  a  12  x  12  grid. 
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Is  not  Intended  to  denigrate  the  extensive  research  efforts  on  Invlscld  PUL 
solvers  and/or  fixed  node  POE  methods;  but  It  does  suggest  that  efforts  to 
accommodate  Laplacian  operators  In  otherwise  Invlscld  solution  methods  and 
efforts  to  Investigate  more  optimal  adaptive  grid  methods  for  use  In  many 
existing  PUE  methods  which  are  applied  to  advectlon-dlffuslon  problems  should 
now  assume  greatly  Increased  significance.) 

It  was  apparent  at  this  stage  of  research  that  our  MFE  results  in  2-0 
were  continuing  a  trend  which  appeared  in  previous  1-0  results.  There,  MFE 
solutions  of  both  the  Navier-Stokes  and  physically  dissipative  continuum 
mechanics  equations  in  1-0  exhibited  perhaps  unprecedented  simultaneous 
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Figure  26.  Isometric  view  of  Burger  s  test  example  at  t  *  1.8  with 
v  *  0.002  on  a  16  x  16  MFE  mesh.  Shock  gradients  have 
magnitudes  of  approximately  103  in  this  solution  for  u. 


resolution  of  extremely  disparate  microscale  and  macroscale  physical  pro¬ 
cesses.  While  2-U  MFE  results  which  emerged  during  this  period  exhibit 

similar  promising  features,  numerous  mathematical  problems  still  require 
resolution  in  order  to  attain  fully  the  desired  levels  of  success  in  truly 
large-scale  physical  problems  in  2-b.  Clues  to  these  problems  can  be  seen  in 
Figure  26.*  For  example,  the  irregularity  of  the  grid  triangles  in  the  face 
of  the  shock  could  eventually  prove  to  be  troublesome.  Similarly,  the  small 
oscillation  at  the  base  of  the  shock  in  this  run  is  unsatisfactory,  even 
though  it  can  be  eliminated  in  any  number  of  ways.  Extensive  testing  and 
analysis  has  indicated  that  the  causal  mechanisms  underlying  such  mesh  irre¬ 
gularities  and  oscillations  in  2-U  can  be  associated  with:  (i)  time  step  and 
error  control  policies  in  the  basic  OUE  integrator  of  Gear  which  is  presently 
used,  (ii)  convergence  properties  of  the  linear  solver,  and  (iii)  limitations 
in  the  first-generation  regularization  functions.  Intensive  investigation  of 
these  factors  reveal  the  needs  in  future  research  which  are  discussed  below. 

UUE  SOLVERS  FOR  PDE  METHODS 

The  current  status  in  this  task  area  is  that  most  existing  OuE  solvers 
are  not  well -suited  for  ready  implementation  in  either  the  MFE  method  or 
numerous  other  advanced  PDE  methods.  This  critical  comment  is,  again,  not 
intended  to  denigrate  the  impressive  advances  in  OUE  research  and  development 
during  the  past  decade;  instead,  it  is  intended  to  bring  a  strong  new  focus 
upon  the  needs  of  PDE  solution  methods,  in  general,  and  more  specifically 
upon  the  pressing  needs  of  adaptive  grid  PDE  methods  which  may  involve  large 
numbers  of  discretized  equations  with  highly  distorted  grids.  Large  distorted 
grid  meshes  may,  in  turn,  augur  for  iterative  linear  solvers  which  can  solve 
poorly  conditioned  matrix  equations,  as  will  be  discussed  further  below. 

A  basic  difficulty  with  most  stiff  OUE  software  has  come  to  the  fore  in 
the  present  MFE  research;  i.e.,  most  ODE  solver  packages  have  been  designed 
to  accommodate  many  different  types  of  classic  OUE  problems.  By  classic,  one 

*It  is  apparent  that  these  suggested  mathematical  problems  will  have  to  be 
resolved  not  only  for  the  MFE  method  but  also  for  most  other  advanced  PUL 
methods  which  may  seek  to  solve  the  difficult  advection-diffusion  equations 
which  frequently  arise  in  physical  problems. 


refers  here  to  such  problems  as  chemical  kinetics  systems  In  which  the  depen¬ 
dent  variables  (e.g.  species  concentrations)  are  all  generlcally  the  same. 

The  error  and  time  step  controlling  policies  In  solvers  for  classic  DDE  pro¬ 
blems  are  usually  less  than  satisfactory  for  applications  of  the  QUt  package 
to  PDE  solution  methods.  In  PDE  systems,  the  spatial  dependence  of  generic- 
ally  dissimilar  variables  comes  into  play.  In  fluid  dynamics  problems,  for 
example,  the  overall  array  of  POE  variables  which  have  been  discretized  on  N 
grid  nodes  (xj,  •••  can  be  represented  as  {p^,  mj,  Ej,  P^.  m^,  Eg,... 

PNmN,  enI»  where  pi  B  p(*i),  p2  E  p(x2^ »  etc.  An  DDE  solver  then  operates  on 
this  array  of  discretized  POE  variables  as  a  single  large  vector  |yit  y^,  y3, 

y4*  •••  >3N-2*  y3N-l*  y3wl*  where  yl  *  pl*  y2  a  ml»  *3  =  El*  *4  s  p2*  etc- 
Because  the  error  control  policies  In  the  Gear  OUE  package,  for  example,  are 

based  upon  L2  norm  of  all  normalized  quantities  / (y ^ )max»  unacceptably  large 

errors  can  be  admitted  in  some  individual  components  of  p,  m,  or  E  at  arbi¬ 
trary  spatial  locations.  A  much  better  measure  for  error  control  policies  in 
PDE  applications  are  maximum  norms  applied  to  each  discretized  PDE  variable. 
The  implementation  of  alternative  norms  is  found  to  extend  deeply  into  the 
logical  structure  of  most  DDE  software  packages,  and  alterations  must  usually 
be  performed  by  someone  who  is  intimately  familiar  with  the  ODE  package. 

In  view  of  such  considerations,  we  devoted  significant  levels  of  effort 
to:  ( 1 )  modifications  of  Gear's  basic  ODE  method  for  MFE  computations;  this 
involved  wholesale  alterations  of  the  Internal  Gear  code  structure  and  also 
extensive  considerations  of  scaling  of  MFE  problem  variables;  and  (ii)  in¬ 
evitably,  the  development  of  entirely  new  ODE  integraton  procedures  which 
better  serve  PDE  solution  needs. 

The  extensive  modifications  to  Gear's  method  have  sufficed  to  solve 
moderately  challenging  PDE's  with  modest  nunbers  of  MFE  grid  nodes  as  was 
seen  above.  But  it  is  now  clear,  also,  that  completely  new  ODE  code  struc¬ 
tures  will  be  needed  in  pending  large-scale  MFE  computations.  We  have, 
therefore,  undertaken  the  development  of  a  low-order  Runge-Kutta  integration 
package  for  MFE  computations.  This  solver  addresses  several  PDE  needs: 

First,  error  control  measures  operate  on  flexibly  ordered  variable  arrays 
using  maximum  norms  on  a  (PDE)  component-by-component  basis.  Second,  PDE 
solutions  have  been  found  to  require  much  more  gradual  time  step  advancement 


policies  than  have  been  built  Into  most  classic  ODE  solvers.*  This  distinc¬ 
tion  between  classic  ODE  and  POE  time  step  properties  apparently  stems  from 
the  fundamental ly  coupled  space-time  dependences  in  POE  systems,  vis  \  vis 
classic  ODE  problems  which  have  no  direct  or  Implied  spatial  dependences. 
Whereas  it  is  computationally  worth  the  effort  to  attempt  very  large  incre¬ 
mental  increases  (sometimes  by  several  orders  of  magnitude)  in  At  in  classic 
ODE  applications  —even  if  such  attempts  may  sometimes  fail— one  finds  that 
the  computational  penalties  for  unsuccessful  large  At  increases  are  much  more 
severe  in  those  POE  applications  where  space-time  couplings  augur  intrinsic¬ 
ally  for  more  gradual  At  advancement  policies.  Third,  time  step  control 
policies  in  the  new  ODE  solver  also  incorporate  convergence  criteria  from 
iterative  linear  systems  solvers.  Such  iterative  solvers  should  henceforth 
be  used  in  largescale  MFE  computations  in  the  interest  of  minimizing  computer 
memory  requirements.  Fourth,  low-order  OUE  methods  are  now  used  because 
high-order  solvers  provide  no  apparent  advantages  in  MFE  applications.  Low- 
order  methods  serve  to  simplify  the  numerical  logic,  improve  the  code  relia¬ 
bility,  and  avert  possible  errors  associated  with  changes  of  order  which  are 
sometimes  present  in  classic  ODE  system  solvers.  Finally,  constraints  on 
allowable  fractional  changes  in  PUE  dependent  variables  are  incorporated  in 
the  overall  time  step  control  policy  in  the  new  ODE  solver. 

LINEAR  SOLVERS  FOR  THE  MFE  METHOD 

Advecti on-diffusion  equations  have  steadfastly  resisted  (if  not  defied) 
satisfactory  numerical  solution  whenever  they  have  been  used  to  describe 
physical  processes  over  highly  disparate  scales.  Such  problems  occur,  for 
example,  in  numerous  applications  of  Navier-Stokes  equations  to  viscous 
compressible  fluids  which  may  contain  shear  layers,  shocks,  and  separated 
flows.  ,  The  basic  difficulty  derives  from  the  nature  of  the  matrix  equations 
which  must  be  solved  in  numerical  PUE  methods  that  are  applied  to  these 
problems.  The  matrix  equations  for  discretized  advection-diffusion  PDE's 
are  large,  sparse  linear  systems  in  which  the  matrices  are  non-symmetric  and 
are  not  necessarily  dominated  by  large  terms  on  the  diagonal.  The  skewness 
of  these  PDE  matrices  can  become  quite  large  for  large  At's  and  for  highly 

*These  policies  also  extend  deeply  into  the  ODE  code  structure. 


distorted  grid  meshes,  both  of  which  are  key  factors  in  efficient  solutions 
of  difficult  advecti on-diffusion  equations. 

The  simple  advecti on-diffusion  equation,  yt  +  £  •  vy  ■  vV2y*  can 
be  used  to  illustrate  some  of  these  features.  Upon  discretization  of  the 
advection-dif fusion  equation,  a  linear  system  of  the  form  (A +3)X  *  C  must 
be  solved.  The  matrix  A  represents  the  stencil  associated  with  the  advection 
operator,  and  the  matrix  B  represents  the  stencil  associated  with  a  nine- 
point  difference  scheme  for  the  Laplacian  operator  in  2-U.  For  represen¬ 
tative  values  of  At,  these  matrices  may  contain  elements  with  the  scaled 
magnitudes  shown  below: 

(0  15  0  \ 

-15  0  15  (3) 

0-15  0  / 

/I  4  1\ 

B  =  I  4  -20  4  J  (4) 

\  1  4  1  /  ,  and 

/  i  w  1\ 

(A+B)  =  -11  -20  19  1  .  (5) 

\  1  -11  1/ 

Testing  and  analysis  has  revealed  that  most  available  linear  solvers  have 
relatively  poor  rates  of  convergence  when  such  significant  large  elements  can 
occur  away  from  the  diagonal  in  non- symmetric  matrices.  Such  iterative  matrix 
solution  methods  as  conjugate  gradient,  multi-grid  and  numerous  other  modern 
linear  systems  solvers  (which  may  work  well  for  symmetric  matrices  in  dis¬ 
cretized  elliptic  equations  and/or  for  uniform  grid  meshes)  do  not  converge 
satisfactorily  in  presently  considered  advection-diffusion  problems.  The 
source  of  difficulty  clearly  derives  from  the  highly  skewed  matrices  and 
dominant  off-diagonal  terms.  We  have  also  shown  that  the  direct  L-U  decom¬ 
position  method  which  has  been  used  in  the  bear  method  until  recently  becomes 
both  noisy  and  computer  storage  limited  when  large  bandwidths  arise  in  pro¬ 
blems  with  more  than  a  moderate  number  of  MFE  grid  nodes. 

Having  identified  more  clearly  the  significance  of  these  issues,  we  have 
developed,  in  collaboration  with  Professor  Keith  Miller,  one  promising  new 


i  r 

i  p 

a 


approach  to  handling  more  effectively  these  imposing  PUE  requirements  on 
linear  systems  solvers.  This  new  matrix  solution  scheme  has,  so  far,  achieved 
good  convergence  rates  for  At's  which  may  be  10  to  20  times  greater  than  the 
large  values  of  t  called  for  by  the  OUE  integrator.  (It  is  generally  hoped 
in  PDE  solutions  that  the  time  step  size  is  determined  by  the  truncation 
error  of  the  OUE  integrator  and  not  oy  severely  limited  convergence  proper¬ 
ties  of  the  linear  solver.)  This  advanced  linear  solver  has  solved  the 
Burger's  equations  discussed  above  with  the  same  CPU  cost  as  the  direct  L-U 
decomposition  method  in  the  Gear  solver,  but  with  greatly  reduced  storage 
requirements. 

REGULARIZATION 


^  Regularization  techniques  have  rarely  been  used  systematically,  if  at 

all,  in  PUE  research  in  the  past.  There  is  thus  presently  great  confusion 
and  misunderstanding  of  the  role  of  regularization  techniques  in  PUE  solu¬ 
tion  methods.  On  the  one  hand,  many  practitioners  of  conventional  PUE  methods 
2  G  suspect  that  regularization  is  an  unfair  trick  by  which  one  can  force  Put 

solutions  to  come  out  in  any  desired  manner  and  that  such  methods  can  there¬ 
fore  not  be  trusted.  On  the  other  hand,  regularization  methods  are  proving 
to  be  valid  and  powerful  mathematical  tools  which  can  now  be  applied  to 
achieve  effective  grid  movement  criteria  systematically  and  to  ensure  that 
high  POE  solution  accuracy  is  also  achieved  in  the  process. 

Only  the  simplest,  first-generation  regularization  functions  have  been 
used  in  2-U  MFE  problems  to  date.  These  penalty  functions  act  like  springs 
and/or  dashpots  in  their  action  on  mesh  triangle  altitudes  (see  Figure  27). 
The  current  penalty  functions  allow  mesh  triangles  to  distort  more  or  less 
arbitrarily,  so  long  as  altitude  magnitudes  remain  positive  and  maintain  some 
tl  designated  minimum  separation.  This  simple  strategy  has  worked  remarkably 

well  in  a  surprisingly  broad  range  of  2-0  problems  considered  to  date,  and  we 
continue  to  probe  the  limits  of  adequacy  of  this  simple  first-generation 
regularization  method.  In  the  Burger's  test  example  discussed  above  with  v  = 
0.002  and  a  16  x  16  MFE  mesh,  minimum  nodal  separations  were  chosen  to  be 
much  smaller  than  v.  As  seen  in  earlier  figures,  the  grid  triangles  can 
become  very  irregular  and  assume  configurations  with  extremely  large  aspect 
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ratios  (0(1Q2)  to  O(103) ) .  Figures  28  and  29  below  show  this  same  Burger's 
test  problem  run  with  "softer"  dashpot  forces  acting  on  triangle  altitudes 
than  In  the  run  shown  previously  In  Figure  26.  The  triangles  in  the  face  of 
the  shock  are  more  regular  In  this  latter  run,  and  their  compaction  in  the 
center  of  the  wavefront  resolves  the  region  where  Intense  shearing  occurs. 
Figure  30  shows  the  projection  of  the  MFE  grid  mesh  on  the  x-y  plane  In  this 
example.  Selected  mesh  connections  have  been  traced  In  heavy  Ink  In  these 
latter  figures  in  order  to  Indicate  the  general  migration  pattern  of  the  MFt 
nodes.  It  is  evident  that  extremely  large  mesh  triangle  aspect  ratios  have 
been  readily  maintained  near  the  shock  front,  still  using  the  improved  L-U 
decomposition  linear  solver  in  the  modified  Gear  integrator  in  both  of  these 
runs.  Perhaps  more  striking  is  the  large  degree  of  fluid  shearing  which  is 
resolved  in  the  face  of  this  skewed  shock  front  by  the  MFt  grid. 

In  the  long  run,  there  is  no  fundamental  reason,  or  desire,  for  the  MFt 
mesh  to  always  be  as  highly  skewed  as  the  physical  flow  lines  in  order  to 
accurately  resolve  such  shear  layers.  It  was  nevertheless  encouraging  at 
this  stage  of  work  that  the  the  MFt  method  is  able  to  handle  such  large  grid 
aspect  ratios  effectively.)  Alternative  regularization  criteria  which  would 
promote  mesh  homogeneity  (e.g.,  by  minimization  of  grid  triangle  aspect 
ratios)  remain  under  Investigation. 

We  have  thus  barely  opened  a  potentially  vast  area  of  investigation  of 
regularization  techniques  for  adaptive  PUE  meshes.  Future  work  should  also 
consider  highly  local  regularization  schemes  which  would  decrease  the  domain 
of  influence  of  current  penalty  functions,  particularly  in  view  of  the  very 
attractive  MFE  properties  which  are  emerging  in:  (i)  alternative  co-ordinate 
systems  and  (11)  the  treatment  of  interface  phenomena. 


Figure  27.  Schematic  Representation  of  First-Generation  MFE 
Regularization  Functions  in  2-0. 


Figure  29.  Isometric  view  of  the  2-D  solution  of  the  velocity 

component,  u,  in  the  Burger-like  example  with  \>  =  0.00 
on  a  16  x  16  MF£  grid.  Note  that  the  viewing  angle  is 
rotated  by  90*  for  a  clearer  view  of  the  doubly  skewed 
wavefront . 
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Figure  30.  MFE  grid  projections  on  the  x-y  plane  in  the  2-0 

Burger-like  example  on  a  16  x  16  grid  with  vs  0.002. 
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ALTERNATIVE  CO-ORDINATE  SYSTEMS 

Initial  work  was  undertaken  on  2-D  MFE  calculations  in  cylindrical  co¬ 
ordinates.  The  results  in  this  area  provide  guidance  for  future  developments 
in  spherical  co-ordinates.  The  apparently  natural  elimination  of  singular¬ 
ities  at  the  origin  by  the  MFE  method  is  of  major  interest  because  such 
singularities  have  been  historically  troublesome  in  many  PDE  methods.  For 
example,  transport  equations  contain  in  cylindrical  co-ordinates  advection 
operators  of  the  form  (ry),  where  r  is  the  radial  co-ordinate;  and  y  is 
a  dependent  variable  of  the  PDE  system.  Singularities  or  other  anomalous 
features  frequently  arise  in  discretized  representations  of  the  term  (y/r) 
as  r— *0.  But  the  MFE  discretization  is  formulated  in  terms  of  well-defined 
inner  prodv^s  which  eliminate  such  possible  singularities.  For  example, 
the  inner  products  of  the  term  (y/r)  with  the  basis  function  o,  taken  over 
the  interval  Ar,  is  given  by 


J(y/r)  •  a  •  rdr  =  Jy  •  adr  . 


The  integral  of  a  *  y  dr  is  essentially  analytic  and  is  readily  eval¬ 
uated  everywhere  on  the  problem  domain.  This  attractive  MFE  property  in 
cylindrical  co-ordinates  obviously  holds  in  a  similar  manner  in  spherical 
co-ordinates.  The  properties  of  these  r-weighted  norms  are  naturally 
different  than  the  MFE  inner  products  which  were  used  in  the  Cartesian 
co-ordinate  systems  considered  in  earlier  MFE  work.  Analysis  and  testing 
of  these  properties  associated  with  r-weighted  norms  and  of  node  controlling 
penalty  functions  in  cylindrical  coordinates  has  been  limited  so  far,  but 
more  extensive  work  should  continue  in  this  area  in  the  future. 


Phase  Transition  Problem:  A  Stefan  Model  of  a  Transiently  Heated  Pipe 


Imbedded  in  Permafrost 


In  a  continuing  analysis  of  MFE  node  movement  properties  during  the 
third  year  of  this  work,  a  basically  cylindrical  problem  was  considered 
using  initially  cartesian  co-ordinates.  This  problem  is  solved  for  the 
transient  melting  of  permafrost  due  to  the  presence  of  a  heated  pipe  buried 
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within  the  permafrost.  The  material  phases  of  Ice,  water,  and  mush  are  des 
crlbed  In  terms  of  their  respective  enthalpies  In  this  Stefan  model.  The 
diagram  below  shows  a  plot  of  enthalpy  versus  temperature  for  the  material 
phases  of  Ice,  water,  and  slush. 

The  variables  and  constants  are: 

h,  enthalpy  per  unit  volume 
f,  temperature 

X,  heat  capacity 
k,  thermal  conductivity 

Y,  latent  heat  of  fusion 

The  POE  for  enthalpy  is 

=  <V2(T(h))  (30) 


(ergs/cm3) 

CK) 

(ergs/ cm3- *K) 
(ergs/cm-sec-*K) 
(ergs/ cm3) 


This  enthalpy  PDE  can  be  rewritten  as 


(31) 


Figure  31.  Diagram  of  Enthalpy  Versus  Temperature  In  Permafrost 


The  generalized  heat  capacity  can  be  obtained  from  the  diagram  above. 

It  Is  given  by 

C(T)  =  |y  =  X+  y  •  6(T)  (32) 

where  <$T  Is  a  sharply  peaked  function.  In  this  example 

S(T)  C(T  -  e){T  ♦  e»2  (33) 

16e5 

Letting  A  =  1,  y  =  2,  and  <  =  1,  the  PDE's  above  become 

[1  +  26(T)  =  V2T  (34) 

3t 

At  t  =  0,  T  =  -1.  everywhere.  Dlrichlet  conditions  are  maintained  on 
the  top  (y  *  1.),  bottom  (y  =  0),  and  vertical  right  boundary  (x  =  1.) 
Symmetry  boundary  conditions  apply  on  the  vertical  left  boundary  (x  =  0.), 
except  at  the  pipe.  A  time-dependent  Dirichlet  condition, 

T  *  -1  +  4  £l  -  ,  (3b) 

describes  the  transient  pipe  temperature  in  the  left  boundary.  That  is,  the 
pipe  is  heated  from  an  initial  temperature  of  -1.0  to  a  maximum  temperature 
3.0.  Figure  32  shows  the  original  MFE  grid,  which  is  essentially  Cartesian. 
After  an  initial  Induction  period,  the  MFE  solutions  exhibit  an  outward, 
circularly  propagating  melt  front  which  is,  in  fact,  an  enthalpy  shock  cor¬ 
responding  to  the  large  gradient  in  enthalpy  at  T  s  0  +  e  associated  with 
the  melting  phase  change.  Figures  33  to  35  show  the  ensuing  evolution  of 
this  enthalpy  shock  using  a  9  x  17  grid  i 

It  Is  Interesting  to  note  in  Figures  33  and  35  that  the  initially  Car¬ 
tesian  MFE  grid  nodes  migrate  readily  to  circular  orientations  before  the 
onset  of  melting.  After  melting  has  begun,  the  MFE  node  contours  track  the 


enthalpy  shockfront  curvature  at  T  *  0  with  high  accuracy.  At  later  times, 
the  MFE  nodes  need  not  be  strictly  conformal  with  the  shock  front,  per  se. 
(See  Figures  34  and  35.)  High  accuracies  are  nevertheless  maintained  due 
to  the  Interpolation  features  and  the  residual  minimization  properties  which 
are  Inherent  In  finite  element  formulations.  Node  biasing  effects  which  are 
associated  with  the  grid  trlangulatlon  are  also  clearly  evident  In  these 
results  at  late  problem  times.  (See,  for  example,  the  right-hand  portions 
of  the  domain;  the  grid  pattern  at  the  top  right  differs  markedly  from  the 
bottom  right  region.  It  is  noteworthy  that  very  little  is  occurring  physic¬ 
ally  in  this  region  which  contains  the  highest  degree  of  grid  biasing.)  It 
is  apparent  that.  In  this  particular  problem,  the  MFE  nodes  exhibit  a  potent 
capacity  to  resolve  fairly  complex  topological  features  with  a  small  number 
of  grid  nodes,  despite  the  presence  of  significant  grid  biasing  features. 

In  the  process  of  obtaining  these  results,  the  following  features  were 
observed  which  suggested  new  MFE  research  in.  Natives : 

(i)  The  numerical  conditioning  of  the  MFE  integrations  can  become 
poorly  conditioned  in  difficult,  highly  complex  flow  regimes.  Numerous 
possible  causes  and  remedies  which  should  be  Investigated  include:  (i) 
improperly  posed  or  incomplete  physical/mathematical  models;  (ii)  highly 
sheared  grid  cells;  ( 1 i  1 )  Implicit  coupling  of  incommensurate  PUE  depen¬ 
dent  variables  which  can  give  rise  to  poorly  conditioned,  non-symmetric 
matrices  which  are  devoid  of  diagonally  dominant  elements;  (iv)  complex 
flow  configurations  near  boundaries;  and  (v)  non-optimal  regularization. 

(ii)  The  simple  grid  triangle  connection  scheme  which  is  used  in  this 
example  (with  all  grid  diagonals  having  the  same  orientation)  can 
exhibit  a  distinctive  biasing  effect  in  the  evolving  MFE  mesh  orien¬ 
tations.  In  many  problems  such  grid  biasing  is  of  no  practical  con¬ 
sequence.  For  example,  MFE  solutions  of  planar  wavefronts  maintain 
extremely  high  accuracies  despite  transverse  grid  biases.)  In  other 
large-gradient  problems,  grid  biasing  has  been  found  to:  (i)  hamper 
computational  efficiency,  and  (11)  possibly  impair  numerical  condition¬ 
ing  and  accuracy.  Several  alternative  grid  connection  schemes  are  logi¬ 
cal  candidates  for  testing  and  Implementation  in  order  to  improve,  and 


possibly  eliminate,  such  grid  biasing  effects.  Alternative  regulariza¬ 
tion  methods  are  also  beginning  to  show  some  promise. 

(ill)  Matrix  solvers.  Extensive  testing  of  matrix  solvers  for  fully 
implicit  MFE  calculations  indicated  three  general  size  distinctions  in 
the  choice  of  matrix  solvers  for  use  on  CRAY-capacity  computers;  these 
are: 

e  Small  problems;  e.g.,  gas  dynamics  equations  —  which  are  discussed 
in  the  next  sample  problem  —  for  the  dependent  variables  p,  u,  v, 
E,  x,  and  y  on  a  10  x  20  grid  mesh,  or  smaller.  Direct  L-U  decom¬ 
position  methods  can  be,  and  have  been,  used  quite  effectively  in 
small  problems.  Other  small  problems  Include  burger's  equations  on 
a  25  x  25  mesh,  as  described  above. 

•  Medium  problems;  e.g.,  gas  dynamics  on  a  30  x  3U  grid  mesh.  Vari¬ 
ous  iterative  multi  grid  methods  exhibit  some  promise  in  this  size 
range.  Additional  research,  development,  and  implementation  is 
very  much  in  order  here. 

•  Urge  problems  are  defined  by  default  from  the  items  immediately 
above.  Practical  solutions  of  large  problems  are  an  ultimate  ob¬ 
jective  for  applications  of  the  MFE  method.  Advanced  iterative 
solvers,  including  perhaps  dynamic  ADI  and/or  other  equation¬ 
splitting  alternatives  remain  to  be  investigated  in  conjunction 
with  both  medium  and  large  problems. 


Gas  Dynamics  Problems  Which  Address  Physical  versus  Non  Physical  (Numerical) 

Dissipation  Effects  in  Gases 

It  is  well  known  in  aerodynamics  that  shockwave-boundary  layer  inter¬ 
actions  are  dependent  upon  local  viscous  dissipation  and  thermal  conduction 
processes  for  the  distribution  of  Internal  and  kinetic  energies.  This  depen¬ 
dence  can  In  many  applications  be  a  sensitive  determinate  of  macroscopic 
flow  properties.  Such  shock-boundary/shear  layer  processes  can  impose  severe 


-44- 


demands  upon  computer  codes  because  the  most  reliable  computations  of  the 
essential  physical  processes  must  first  Include  the  physical  dissipation 
operators  of  the  full  havler-Stokes  equations  and  then  use  optimally  dis¬ 
tributed  grid  densities  which  can  accurately  resolve  the  actual  physical 
dissipation  processes  and  eliminate  the  numerical  dissipation  processes. 

To  gain  a  full  appreciation  of  the  origin  and  nature  of  physical  dlsslptlon 
processes.  It  Is  well  to  recall  the  development  of  the  Kavier-Stokes  equa¬ 
tions  In  the  context  of  1-D  kinetic  theory.  That  Is,  the  fluid  equations 
can  be  written  as: 

&  *  h (p,)  * 0  iM 


3(pv)  4.  3  ,2)  *  3  n(0)  +  nU)  +  J2)  + 

Tt~  3x  { pv  J  Tx  p  +p  +p  +*“ 


3E  ,  3  /r  .  3  ,  ,  3 

st  *  s7  ,E,)  '  "57  (pv)  0 


(o)  .  (i)  .  u)  . 

+  q  +  q  +  . . . 


A  zero-th  approximation  of  the  kinetic  theory  for  gases  uses  the  consti¬ 
tutive  relations  for  an  inviscid,  non-conducting  fluid;  i.e., 

(0)  pRT  ,  iwr  1  2. 

p  ■  p  =  ■  (y  -  1)  (E  -  2PV  )  (39) 


,  ■  q(0)  =  0 


These  relationships  yield  the  well-known,  inviscid  Euler  equations. 

A  first  approximation  of  the  kinetic  theory  for  gases  gives  the  Navier- 
Stokes  equations  according  to: 

p  =  p(0)  +  p(1)  «  (Y  -  DU  -  ±pv2)  +  (41) 


"  ■  «<0)  * 1,(11  -  - 


A  second  approximation  of  kinetic  theory  gives  the  Burnett  equations,  etc. 
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In  order  to  address  the  primary  Issue  of  the  potential  significance  of 
physical  dissipation  processes  as  determinates  of  transient  macroscopic  flow 
properties,  we  wish  to  distinguish  the  relative  roles  in  code  calculations 
of  numerical  vis  a  vis  the  physical  dissipation  effects  which  appear  on  the 
right  hand  sides  of  Equations  (37}  and  (38).  For  this,  we  have  solved  a 
simple  shock-boundary  layer  reflection  process  in  1-U.  This  example  is 
frequently  referred  to  as  the  anomalous  wall  heating  problem. ^  In  this 
problem,  anomalously  high  temperatures  are  calculated  by  many  existing  shock 
hydrodynamics  codes  for  the  reflection  of  a  planar  shock  from  an  infinitely 
reflecting  wall  in  slab  geometry.  The  MFE  results  which  follow  indicate  that 
the  anomalous  aspects  are  eliminated  when  numerical  dissipation  is  suppressed 
and,  most  importantly,  when  the  physical  dissipation  processes  in  the  Navier- 
Stokes  equations  are  accurately  resolved  in  the  transient  reflection  process. 
The  following  test  problem  in  1-0  slab  geometry  illustrates  these  results; 


Initial  Conditions: 
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0.  _<  x  _<  2. 
0.  <  x  <  2.0 
Ax0  <  x  <  2.0 
0  ?  x  7  Ax0 
x  =  0. 


Boundary  Conditions: 

Reflection  at  x  *  0. 
Dirichlet  at  initial  values 


at  x  «=  2.0 


Rankine-Hugoniot  Solutions  for  Infinite  Shock  (t-*»): 


s  =  1/3 
d+  =  4.0 
i*  =  0.5 
*+  =  0. 

3+  =  1.33 


p-  =  1.0 

e“  *  0. 
v-  =  -1. 

p~  e  0. 


The  time  evolution  of  this  shock  was  solved  by  the  MFE  method  in  two 
ways:  First  the  full  Navier-Stokes  equations  were  solved  accurately  using 
alternative  values  of  v  ■  4y/3  and  *c.  These  solutions  are  denoted  by  N-S  in 
the  accompanying  figures.  In  one  set  of  N-S  solutions,  v  *  k  *  0.01,  which 
is  unrealistically  large  but  which  permits  comparisons  to  other  fixed  node 


POE  solutions  that  may  use  on  the  order  of  100  to  several  hundred  grid  nodes. 
In  another  set  of  N-S  solutions,  v  *  <  *  0.001,  which  approximates  physically 
realistic  values  for  actual  dissipation  processes  in  gases.  Second,  the 
variable  denotes  internal  energy  per  unit  mass  in  the  accompanying  figures, 
and  anomalous  diffusion  (denoted  by  A  in  the  accompanying  figures)  was  simu¬ 
lated  by  including  a  diffusion  term  vrpxx  in  Equation  (3b).  This  effectively 
simulates  some  form  of  uncontrolled  numerical  diffusion  which  is  present 
intrinsically  in  all  numerical  POE  solution  methods.  (Uncontrolled  numerical 
diffusion  is  the  only  source  of  dissipation  in  numerical  solutions  of  the  in- 
viscid  Euler  equations.)  In  the  results  which  follow,  we  have  verified  that 
the  MFE  solutions  of  the  Navier-Stokes  equations  have  reduced  all  numerical 
or  other  anomolous  diffusion  effects  to  imperceptible  low  levels  and  that  the 
observed  behavior  of  shock  interactions  is  associated  with  physical  dissipa¬ 
tion  operators. 

At  t  *  U+,  the  shock  incident  on  the  origin  is  in  the  incipient  state  of 
outward  reflection.  At  t  =  0.0b,  Figures  36  and  37  show  that  the  calculations 
of  the  reflected  shock  with  uncontrolled  diffusion  (or  simulated  numerical 
diffusion)  tend  immediately  to  overheat  in  e  and  to  correspondingly  under¬ 
shoot  in  p  relative  to  the  Navier-Stokes  solutions.  Although  these  transient 
solutions  are  not  near  their  steady  state  values  at  this  early  time,  it  will 
be  seen  that  the  ensuing  evolution  toward  equilibrium  is  quite  sensitive  to 
both  the  magnitudes  and  the  nature  of  the  dissipation  processes  in  the  compu¬ 
tations. 

Figure  38  shows  that  at  t  *  0.15,  the  lip  of  the  shock  in  the  Navier- 
Stokes  solution  is  approaching  the  steady  state  value  of  P  «=  4.0,  and  the 
anomolous  dissipation  solution  lags  by  a  significant  margin.  The  fluid 
buildup  at  the  front  of  the  shock  is  evident  here  because  the  fluid  near 
the  origin  has  stagnated  while  additional  fluid  continues  to  stream  in  toward 
the  origin  from  the  region  to  the  right  of  the  shock.  Figure  39  shows  that 
the  anomolous  dissipation  results  continue  to  lag  behind  the  Navier-btokes 
solutions  to  a  significant  degree  at  t  s  0.300.  At  t  *  2.0,  the  Navier-Stokes 
solutions  have  approached  steady  state  Rankine-Hugoniot  conditions  (not  shown 
in  Figure  40),  and  the  anomalous  dissipation  solution  has  still  not  reached 
the  Rankine-Hugoniot  values  in  the  vicinity  of  the  origin.  The  anomalous 


Figure  38.  MFE  solutions  of  the  Navier-Stokes  Figure  39.  MFE  solutions  of  the  Navier-Stokes 

equations  and  anomalously  dissipative  equations  and  anomalously  dissipative 

equations  at  time  =  0.150.  equations  at  time  =  0.300. 


wall  heating  effects  due  to  uncontrolled  dissipation  In  the  density  equation 
have  thus  persisted  to  very  long  times  vis  a  vis  the  accurate  solutions  of 
the  Navier-Stokes  equations.  Non-physical  dissipation  which  is  sometimes 
Introduced  as  artificial  viscosity  terms  in  fluid  calculations  can  be  shown 
to  have  similar  anomalous  effects,  as  can  the  nunerical  dissipation  and/or 
numerical  uncertainties  which  are  present  in  solvers  of  inviscid  Euler 
equations. 

Figures  41  and  42  present  the  results  of  another  test  of  sensitivity  of 
the  Navier-Stokes  solutions  to  non-optimal  grid  locations.  In  this  test 
problem,  a  physically  realistic  value  of  v  =  U.Q001  is  used  in  MFE  solutions 
of  the  Navier-Stokes  equations.  We  have,  however,  deliberately  constrained 
the  MFE  grid  nodes  in  this  test  case  so  that  they  do  not  migrate  to  their 
truly  optimal  locations,  as  in  the  results  considered  previously.  Figure  41 
shows  several  significant  features:  (i)  the  shock  gradients  associated  with 
v  *  0.0001  are  extremely  large;  the  accurate  resolution  of  these  gradients 
would  require  several  thousand  nodes  if  a  fixed  node  POE  solution  method  were 
to  be  used,  (ii)  the  Rankine-Hugoniot  solutions  are  approached  much  more 
rapidly  for  the  physically  realistic  values  of  v  than  for  the  larger  values 
of  v  which  are  typically  used  either  tacitly  or  explicitly  in  many  other  POE 
solution  methods,  and  (iii)  the  slight  constraint  on  node  movements  and  thus 
on  nodal  positions  do  not  show  up  immediately;  but  once  the  perturbation 
becomes  significant  (as  seen  in  Figure  42),  its  effects  can  grow  rapidly.  In 
summary,  these  results  demonstrate  that  reflected  shock  solutions  can  be  very 
sensitive  to  non-physical  dissipation  effects  and  to  slight  deviations  from 
optimal  grid  node  positioning,  even  in  adaptive  gridding  methods.  All  of  the 
results  in  this  section  were  obtained  with  approximately  30  MFE  nodes.  As 
many  as  61  MFE  nodes  were  used  to  verify  that  the  MFE  solutions  were  in  fact 
converged  solutions.  We  have  also  compared  the  steady  state  MFE  solutions  of 
fluid  velocities  with  the  analytic  Navier-Stokes  solutions  of  the  steady 
shock.  Excellent  agreement  was  obtained  between  the  MFE  solution  and  the 
analytic  Navler-Stoke  solution.  In  contrast,  it  was  found  after  extensive 
attempts  that  computed  solutions  which  contained  nominal  degrees  of  numerical 
or  other  anomolous  diffusion  could  not  be  forced  by  parameter  manipulations 
to  agree  with  the  analytical  shock  solutions.  In  still  other  MFE  calcula¬ 
tions,  we  have  computed  in  1-0  the  interaction  of  a  shock  with  an  internal 


shear  layer  associated  with  a  contact  discontinuity  ~  In  direct  analogy  to 
the  Irregular  shock  reflection  mechanisms  which  occur  In  higher  dimensions. 
Here  also,  the  computed  macroscopic  flow  properties  were  sensitive  to  accurate 
resolution  of  the  physical  dissipation  processes  In  the  full  Navier-Stokes 
equations.  In  the  absence  of  the  stringent  tests  of  convergence  which  were 
applied  here,  it  can  be  extremely  difficult  to  discern  physical  oscillations 
and  dissipation  effects  from  non-physical  and/or  purely  numerical  oscillations 
and  dissipation  effects.  We  have.  In  fact,  found  It  generally  impossible  to 
simulate  the  effects  of  the  actual  physical  dissipation  processes  in  shock¬ 
boundary  layer  interactions  by  the  use  of  either  artificial  or  parametrically 
controlled  numerical  diffusion  processes  in  Euler  equation  solutions. 

Shjck-on-Wedge  Problem  in  2-D 

A  shock-on-wedge  problem  which  can  be  verified  against  the  experimental 
data  of  Glass  and  co-workers(lu)  has  been  addressed  in  2-D  MFE  calculations. 
This  problem  is  posed  in  rectangular  coordinates,  as  shown  by  the  schematic 
view  of  an  initial  grid  mesh  configuration  in  Figure  43. 


x 

Figure  43.  MFE  Initial  Grid  Zone  Configurations  for  Shock-on-Wedge  Problems. 

In  2-0,  we  have  used  the  mass,  momentum,  total  energy  representation  (p, 
m,  E)  of  the  fluid  equations.  The  Navier-Stokes  equations  in  the  (p,  m,  E) 
representation  are  given  by 


+  £  m 
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■^  +  7  *  (m  m/p)  *  4,  *  (44) 

3E  V' 

~  +  V  •  (Em/p)  -  -£  *  Q  +  V  *  (2-.  *  (m/p))  ,  (4b) 

a  t 

where ^  is  the  generalized  stress  tensor  and  £  is  the  heat  flux  vector.  The 
heat  flux  vector  is  frequently  given  by 

Q  =  -  k£  T  ,  (4b) 

where  T  is  the  temperature.  The  general  stress  tensor  is  frequently  given  by 

=  -pl^  +  t  ,  (47) 

where  I  is  the  identity  matrix.  For  an  ideal  gas  (a  Newtonian  fluid)  in 
Cartesian  co-ordinates 


p  =  p^  =  (Y  -  1) (E  -  m^/2P),  and 


3(m/p)x  2 

Txx  =  2W  ST”'  3y<I  * 

3(m/p)y  i 

xyy  *  2y  3y  -  3  V  (V  •  (m/p)) 


3(m/p)x  3(m/p)y 

Txy  ~  y  3y  +  3x 


It  is  important  to  note  in  Equation  (44)  that  the  quantity  mm  is  a 
dyadic.  Special  attention  has  been  devoted  to  exacting  evaluations  of  those 
inner  product  terms  which  involve  factors  of  -jj  for  small  p.  These  factors 
are  present  in  all  standard  representations  of  the  fluid  equations. 


KV 


ri 


This  2-0  test  problem  considers  the  case  of  regular  shock  reflection  in 
argon  at  Mach  number  =  2.05,  p0  =  150  torr,  p0  •  3.23  x  10“4  g/cm3,  T0  « 
297. 6*K,  and  wedge  angle  0W  =  60*  (see  Figure  44). 
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Figure  44.  Schematic  representation  of  regular  reflection  in 
Shock-on-wedge  experiments. 
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Initial  conditions  calculated  from  the  Rankine-Hugoniot  relationships  are: 


p+  = 

3.23  x 

10"4  g/cm3 

m+  = 

0 

E+  = 

[p+  e+ 

+  m+/2p+]  = 

P-  = 

7.87  x 

10“4  g/cm3 

m_  « 

30.559 

gm  -  cm/sec 

E_  = 

[p_  e- 

+  m£/2p-]  = 

3.0  x  105  ergs/ cm3  and 


2.132  x  106  ergs/ cm3 


Boundary  Conditions 

A  significant  point  of  fundamental  physics  became  manifest  in  the  study 
of  boundary  conditions  for  shock -on-wedge  calculations.  For  a  variety  of 
reasons,  conventional  shock  codes  tend  to  use  slip-type  boundary  conditions. 
In  actuality,  the  velocity  of  air  molecules  is  zero  at  the  types  of  surfaces 


under  present  consideration.  For  the  sake  of  common  comparisons,  and  to  test 
the  sensitivity  of  MFE  shock-on-wedge  computations  to  alternative  boundary 
conditions,  we  first  attempted  to  solve  the  regular  shock-on-wedge  test  pro¬ 
blem  with  slip-type  boundary  conditions,  in  which  the  normal  component  of 
velocity  along  the  bottom  horizontal  surface  and  inclined  wedge  are  zero  and 
the  tangential  velocity  components  are  treated  by  zero-Neumann  (or  symmetry) 
conditions  (see  Figure  44).  These  symmetry  boundary  conditions  are  applied 
similarly  to  the  MFE  grid  nodes  so  that  the  grid  nodes  would  slide  freely 
along  both  the  horizontal  surface  and  the  inclined  wedge  surface.  But  the 
NFE  method  has  always  been  found  to  be  extremely  sensitive  in  its  accuracy 
and  consistency  requirements  to  the  presence  of  any  physical  or  mathematic¬ 
ally  ill-posed  condition.  Here,  the  surface  normal  is  not  defined  uniquely 
at  the  front  of  the  wedge.  Consequently,  the  basic  consistency  requirements 
in  the  MFE  method  impeded  the  numerical  integration  process  because  the 
symmetry  boundary  condition  at  the  foot  of  the  wedge  is  ill-posed,  and  thus 
could  not  be  properly  resolved  as  a  truly  rigorous  PUE  solution.  The  MFE 
method  also  tends  to  be  less  tolerant  than  most  PUE  methods  of  computational 
swindles  which  are  frequently  used  to  accomplish  slip-type  boundary  condi¬ 
tions.  It  thus  became  quite  apparent  that  it  is  possible  to  enforce  slip- 
type  boundary  conditions  only  by  accepting  erroneous  numerical  solutions  in 
the  boundary  layer  region  near  the  front  corner  of  the  wedge.  It  was  further 
apparent  that  these  errors  could  propagate  to  regions  well  away  from  the 
local  source  of  difficulty.  In  view  of  these  results,  we  next  attempted  to 
take  a  large  step  forward  in  terms  of  both  the  physics  and  computations  of 
the  physically  required  non-slip  boundary  conditions.  This  turned  out  to  be 
a  large  step  because  this  shock-on-wedge  problem  has  a  turbulent  boundary 
layer  in  the  duct  represented  by  our  problem  domain.  Nevertheless,  some  key 
points  became  manifest  In  this  task  which  will  be  described  Immediately  below. 

The  initial  shock  conditions  are  shown  in  Figure  45.  This  problem  is 
solved  in  a  1  cm  duct.  For  viscosity  y  «  5  x  ICT^,  the  Reynolds  number  is 
approximately  5.7  x  1Q4.  Such  a  flow  is  clearly  turbulent.  The  MFE  solutions 
of  the  laminar  Navler-Stokes  equations  were  nevertheless  attempted  as  a 
numerical  experiment.  These  results  appear  in  Figures  46  to  55  with  self- 
explanatory  captions.  At  t  *  U.30  the  shock  is  just  starting  to  encounter 
the  wedge  front;  and  the  Isodensity  contours  early  in  the  reflection  process 
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U  *  6.58  x  104  cm/ sec 


P+  *  3.23  x  10"4  g/cm3 
p+  «  150  torr 
v+  =  0. 

E+  *  3.00  x  10®  ergs/cm 


C+  *  3.21  x  104  cm/sec 


P_  =  7.54  x  10“4  g/cm3 
p_  =  760  torr 
v  *  3.76  x  10^  cm/sec 
E  =  2.03  x  10®  ergs/cm3 


Figure  45.  Initial  Conditions  for  Regular  Reflection  of  Planar  . 

Shock  in  the  Experiments  of  Deschambault  and  Glass/10' 


Figure  46.  MFE  Solutions  of  Laminar  Navler-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  a  2.05,  e  =  60°,  p+  =  150  torr,  p+  =  3.23  x 
10’4g/cm3  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 


Figure  47. 


t  *  0.30 


Figure  48. 


MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  (Ms  =  2.05,  0  =  60°,  p+  =  150  torr,  p+  =  3.23  x 
10-4g/cm3  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 


MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  =  2.05,  9  =  60°,  p+  =  150  torr,  p+  =  3.23  x 
10'Vcm^  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 


Figure  52. 


MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  [Ms  =  2.05,  0  ■  60°,  p+  =  150  torr,  p+  =  3.23  x 
10"^g/cnw  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 


Figure  54. 


MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  (M5  *  2.05,  0  =  60°,  p+  =  150  torr,  p+  =  3.23  x 
10_4g/cm3  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 
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Figure  55.  MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  a 
Regular  Shock  Reflection  Experiment  of  Deschambault  and 
Glass.  (M5  =  3.06,  0  =  60°,  p+  =  150  torr,  p+  =  3.23  x 
10-4g/cmJ  in  Argon.)  The  MFE  Grid  is  9  x  51  Nodes  and 
Non-Slip  Boundary  Conditions  are  Used. 
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Figure  56.  MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  an 

Experimental  Plane  Shock  of  Deschambault  and  Glass  Reflecting 
Against  a  Vertical  Wall.  (Ms  =  2.05,  0  =  90°,  p+  =  150  torr, 
p+  =  3.23  x  10-4g/cm3  in  Argon.)  A  3  x  31  MFE  Grid  is  Used. 
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Figure  59. 
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Figure  60. 


MFE  Solutions  of  Laminar  Navier-Stokes  Equations  for  An 
Experimental  Plane  Shock  of  Deschambault  and  Glass  Reflecting 
Against  a  Vertical  Wall.  (Ms  =  2.05,  6  *  90  ,  p+  =  150  torr, 
p+  =  3.23  x  10"4g/cm3  in  Argon.)  A  3  x  31  MFE  Grid  is  Used. 


at  t  *  0.35  have  a  regular  profile.  At  t  *  0.44  the  formation  and  subsequent 
shedding  of  eddies  are  becoming  apparent.  These  are  formed  and  evolved 
according  to  the  physically  inadequate  laminar  dissipation  processes  in  the 
present  Navier-Stokes  solutions  and  not  by  numerical  dissipation.  Analysis 
of  vorticity  generation  rates  indicates  clearly  that  a  turbulent  model  of  the 
physical  dissipation  in  the  boundary  layer  is  required  in  order  to  model  this 
system  according  to  physical  reality.  We  nevertheless  seem  to  be  computing 
the  early  onset  of  turbulent  behavior. 

Mostly,  this  example  shows  an  apparent  capacity  of  the  MF£  method  to 
attempt  to  resolve  extremely  micro-scale  dissipation  processes  simultaneously 
with  the  macroscopic  flow  features.  In  this  sample  problem,  it  was  undoubt¬ 
edly  the  physics  model  which  requires  improvement  more  than  the  numerics. 
Because  local  turbulent  dissipation  rates  exceed  laminar  rates  by  large  fac¬ 
tors,  it  is  likely  that  this  experimental  test  example  imposes  more  severe 
demands  upon  the  PUE  method  than  will  be  encountered  in  more  realistic  physi¬ 
cal  models  of  real  shock  environments  (which  have  much  larger-than-laminar 
dissipation  rates). 

In  order  to  verify  that  the  complex  behavior  in  this  example  was  truly 
associated  with  vorticity  generation  in  the  boundary  layer,  this  same  inci¬ 
dent  shock  was  propagated  into,  and  reflected  from,  a  vertical  wall.  Fluid 
shears  do  not  occur  in  this  example,  and  the  MFE  solutions  appear  in  Figures 
56  to  62.  The  shock  profiles  remain  perfectly  planar  (to  many  significant 
figures)  because  only  compressive  and  expansive  forces  act  in  this  reflection 
process,  neither  numerical  dissipation  nor  anomalous  vorticity  effects  were 
perceptible  in  these  MFE  solutions.  It  was  also  validated  that  the  Rankine- 
Hugoniot  and  Taylor  solutions  of  the  steady  reflected  shock  profiles  were 
obtained  in  these  MFE  solutions. 

SUMMARY 

The  report  above  describes  all  significant  results  of  this  three-year 
investigation  of  the  basic  MFE  method  and  its  properties,  both  the  promise 
and  the  problems  which  may  require  extensive  additional  research  have  been 
indicated. 
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